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Abstract 

The effects of plasma ion motion in self-modulated plasma based accelerators is examined. An 
analytical model describing ion motion in the narrow beam limit is developed, and confirmed 
through multi-dimensional particle-in-cell simulations. It is shown that the ion motion can lead to 
the early saturation of the self-modulation instability, and to the suppression of the accelerating 
gradients. This can reduce the total energy that can be transformed into kinetic energy of accel- 
erated particles. For the parameters of future proton-driven plasma accelerator experiments, the 
ion dynamics can have a strong impact. Possible methods to mitigate the effects of the ion motion 
in future experiments are demonstrated. 
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Plasma based accelerators (PBA) [JJ sustain large amplitude waves that can trap and 
accelerate particles to high energies in distances more than three orders of magnitude shorter 
than conventional devices [2J. Currently, PBAs use < 10 J laser beams [3 J (LWFA), and 
< 1 kJ electron and positron bunches [5H7] (PWFA) to accelerate ~ 1 GeV (LWFA) to 
~ 100 GeV (PWFA) electron bunches. The use of 100 kJ proton bunch drivers to generate 
TeV-class electrons was recently proposed in the so called proton-driven plasma-wakefield 
accelerator (PDPWFA) [8] which relies on very short proton bunches [o z < 10 /im), with 
more than 10 11 protons per bunch to excite plasma waves in the blowout regime [H]. 

The shortest proton bunches currently available, however, are long, with o z > 10 cm. 
Thus, a future PDPWFA experiment pU] will operate in weakly relativistic regimes where 
two stream-like (self-modulational) instabilities [UJ dominate the proton beam dynamics 
and wakefield excitation [TJJ [13]. There are many analogies between self-modulation (S- 
M)of lasers and particle beams [14J, and this scenario is also of relevance for the propagation 
of intense plasma streams in astrophysics [To] . 

Recent work on the S-M of long particle bunches JTUJ [TJJ revealed that the wake am- 
plitude and phase velocity can be nearly constant once the instability saturates. How- 
ever, it was assumed that the background plasma ions are immobile. This will certainly 
not be the case for very long beams when u P i<j z /c ^> 1 (oo P i = a/ 'Annoe 2 /rrii is the ion 
plasma frequency, e the elementary charge, and m 8 the ion mass). It is therefore im- 
portant to understand the role of the ion motion in the S-M of very long beams, i.e. 
Upi&z/c — 1-88 x 10 2 (cr z /10 cm) ^m e /mi^n [cm -3 ] , where m e is the electron mass, and 
to determine conditions for which the plasma ion motion can be minimized. 

Previous work addressed the role of ion motion due to the space charge forces of ultra- 
intense driving electron bunches in the PWFA nonlinear blowout regime [18]. Here we 
explore the motion of the plasma ions in the linear regime, where the ions respond to the 
plasma wave ponderomotive force. 

In this Letter, we show that the background plasma ion motion can play a key role in the 
S-M of proton beams. An analytical model, confirmed by multi-dimensional particle-in-cell 
(PIC) simulations in OSIRIS [17], demonstrates that ions move due to the nonlinear plasma 
wave ponderomotive force on plasma electrons. Future PDPWFA experiments will operate 
in the narrow beam limit for which there has been little work, and where fluid theory may fail 
due to phase mixing. We thus used the plasma sheet model [TB] to describe wake excitation 
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in this limit and derive exact expressions for the ponderomotive force. Full-scale simulations 
in conditions of the PDPWFA reveal that the ion motion suppresses the transverse self- 
modulation instability (SMI), and reduces the accelerating wakefields throughout most of 
the proton bunch. Preventing the ion motion is therefore crucial, in particular for a future 
PDPWFA experiment, and heavy ion plasmas (e.g. Ar + ) are required to prevent it. 

We start with the plasma fluid equations. The linearized continuity equation for the 
ions is dun/ dt + n V ■ v. a = 0, where rij = n + rin is the ion density, n the background 
plasma density, rin <C n and v ;i <C c the perturbed ion density and velocity respectively, 
and t the time. We assume cold ions (Tj = 0) such that the linearized Euler's equation is 
dvji/dt = ZeE/rrii, where E is the electric field, and Z the ion charge. To obtain E we 
ignore the electron inertia and assume that n» = n e which is true if X 2 ,V 2 (n e /rio) <C 1, giving 
eE = /cT e V(nj/-n,o) + F p where F p is the nonlinear ponderomotive force on the electron, T e 
the electron temperature, = \/kT e j 'Anrioe 2 the Debye length, k the Boltzmann's constant, 
and V 2 the Laplacian operator. Differentiating the ion continuity equation, substituting the 
dv^/dt from the Euler's equation, and using E from the electrons' Euler equation then 
gives: 

n a = -n ZV ■ F p , (1) 

where c s = a/ ZkT e /rrii is the ion sound speed, and where (£ = z — ct, r = r, r = t). 
Equation (II]) describes the evolution of the ion density in a warm plasma with rin "C no- 
For a plasma wave with frequency u p , then F p = e 2 /(4m e w 2 )VE 2 |19j . where E is the 
envelope of the plasma wave electric field. 

We consider cylindrically symmetric drivers such that E, n e , and rii only depend on 
the radius r, on z, and on t. In general both radial and longitudinal components of F p 
should be included J2JJ. For cigar-shaped drivers with o z 3> a r , however, |F p j_| 3> |F P ||| 
since d z ~ d r ~ In addition, for narrow beams with a r c/u p , E z E r since 

the motion of the plasma electrons is preferentially along r and not along z. Under these 
assumptions, and for a cold plasma, Eq. ([T]) reduces to: 

2^ n Ze 2 2p2 

Furthermore for narrow and small amplitude wakes the plasma response is dominated by 
the electrostatic forces for E r , so we can obtain F p for narrow bunches using Dawson's sheet 
model [16] where non-relativistic, electrostatic radial oscillations of cylindrically symmetric 
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electron rings are considered. In the absence of trajectory crossing the equation of motion 
for an electron ring is |16j : 

dV _ _c^r u 2 p r 2 e E* 

d£ 2 2 2r m e ' U 

where r is the radial position of the electron ring, and r its initial radial position. The 

first term on the right hand side of Eq. ^ corresponds to the ion channel focusing field 

(El), the second to the electrostatic repulsion from the plasma electrons (E*) inside the ring 

(the charge inside the ring of electrons remains fixed if the electron rings do not cross), and 

the third to the radial force associated with a cylindrically symmetric particle beam with 

electric field E^. We let r = r + Ar(r ,£) and expand the right hand side of Eq. ^ into 

powers of Ar/r yielding: 

2 d 2 Ar _ eE r b (r ,Q _ 
d£ 2 " m e 1 ' 

where the derivatives are evaluated at r = r$. Equation Q is also valid in 2D slab geometry. 
We consider a flat top driver in £ with length a z . Inside the beam (0 < £ < a z ), Ar = 
A< (cos0 - 1) where A<(r ) = ^/f^i^o) > and = {u p £/c)^l + (e/w 2 m e )V r0 £ r b (r ). 

1 /2 

Behind the driver (£ > a z ), Ar = v4>cos0, where v4>(r ) = [Ar 2 + c 2 /u; 2 (dAr/d£) 2 ] , 
Ar = r(cr 2 ) — r , <t> = uj p ^/c + (fio, and (fi is the phase at £ = a z . These trajectories are in 
agreement with PIC simulation results using narrow drivers (cr r <C c/u p ). 

To determine E r (r) self-consistently the expression for the electron trajectories is in- 
serted in E r (r (r,^),r) = E % r (r (r, £),r) + E* (r (r, £), r) + £y(r) if r (r, £) was known. 
To determine r (r, £), r = r + Ar(r ,£) is inverted using the Lagrange's implicit func- 
tion theorem. To lowest order we find r (r, £) by Taylor expanding r(r ) in powers 
of Ar yielding r ~ [r — y4<(cos0 — 1)] — [1 + V r A<(cos0 — l)]Ar for £ < <j 2 , and 
r ~ r — v4> cos + (1 + V r A> cos 0) Ar for £ > cr^, and then by solving these expressions for 
tq. Inserting the result into the right hand side of Eq. Q yields: 

Ef<" = ^<(-l + cos0) (5) 

1 + V r eE r< /(m e u:l) (cos — 1) 



for £ < cr z and : 



E r > cos <t> , \ 

E^ = ^ " , (6) 

1 + V r e-E r> / (m e o; 2 ) cos 
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for £ > a z , where E r ^ = m e u 2 A^/e is the amplitude of the radial plasma wave. When 
E r ^(e/m e u p ) <C 1, and V ' r E r ^(e/m e ujp) <C uj p /c, E r is purely sinusoidal. However, when 
narrow drivers are used, V r E r ^(e/m e ujp) > w p /c and the wake becomes anharmonic even 
though E r ^(e/m e Up) <C 1. Excellent agreement was found between the analytical model 
and PIC simulations [22] using narrow drivers, where linear fluid theory fails to accurately 
describe wake excitation. 

The average E r over a plasma oscillation is finite. This can be interpreted as the 
electrostatic force from the ions that is required to balance a nonlinear average force 



pushing electrons outward. Thus, F t 



p±\ 



is simply obtained by averaging E r over 



one oscillation, i.e. e{E r ) 



Using the fact that J„ 27r [1 + x(cos0 — 1)] 1 d<fi 



-(1 + 2s)-V>, £ 
Jq W cos 0/(l + a; cos < 



cos< 



[1 + ,x(cos0 - 1)] 1 d(f) = (1 + x - yT+ 2x)/(xy/l + 25), and 
.(f) — [1 — (1 — x 2 )~ l l 2 \/x [20], then the average value of Eq. (5) for 



£ < a z , and near the axis becomes : 



{Er) 



m e UpT 



2e 



■V r E, 



2 V r ±J r < 



+ E b r 



(7) 



and the average value of Eq. (|6|) for £ > a z at any r becomes 



(E r ) 



E, 



r> 



V r E. 



r J -'r> 



-.{VrE,, 



(8) 



r>J 



Equations (7l8) generalize the fluid theory nonlinear ponderomotive force (F p = e(E r )) 
for narrow plasma waves and are correct to lowest order in r — r , i.e. they are valid as 
long as E r ^ is calculated for sinusoidal trajectories. The corrections due to anharmonic 
oscillations in Ar could also be included by considering higher order terms in Eq. (Q. To 
make a connection to the usual expression for F p we make a Taylor series expansion for small 
V r E r in Eq. (8) which yields F p± = -e 2 /{Am e ul)V r E 2 + O (PA. Equations (7 8) differ 
from Fp ± when the wakefield becomes nonlinear; the wake can be anharmonic even when 
e/(m e u p c)E r <§C 1 if e/(m e u^)V r E r < 1, if a r <C c/u p [22J. This is the limit of relevance for 
the PDWPFA. 

The ion density perturbations in the wake of proton bunches can now be found by 



substituting (E r ) into Eq. (2). For £ <C c/u P i this yields rii = 
Near the axis where E r ~ rV r .Ey|o, then, for £ < o z 



n 



1 + HfiV ■ (Er) 



n iQ [l 



£.<&z 
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- (1 - -^V r E r <)- 1/2 + ^V r E b r )\. For f > a z , at any r, nf ><Tz 



n 



2 C 2 



Thus, F p x pulls the plasma ions to- 



i - (i - (i - e 2 (v r 4>) 2 MH 4 r 1/2 

wards the axis. Ion density voids are also formed near the edge of the wakefield. 

These theoretical predictions were compared with 2D slab geometry and 3D PIC 
simulations. The simulations used a 3D (2D) window moving at c, with dimensions 
200x4x4(c/u;p) 3 (400x4(c/u;p) 2 ), and a grid with 2560x512x512 (10240x1024) cells in the 
longitudinal and transverse directions respectively. Each cell contains 2x1x1 (2x2) elec- 
trons and ions with mass mi = 1836 m e . The simulations considered = a b rexp [— r 2 /er 2 ]. 
Both finite flat-top and infinite drivers were used. The values for a h and oy are given in Fig. [TJ 
Figure [l] shows the simulations and the predictions from the analytical model (Eqs. Q, ( |7f8| )) 
(where V 2 = d 2 /dr 2 was used in 2D). There is excellent agreement for this range of param- 
eters. Additional simulations with smaller oy confirmed that fluid theory underestimates 
the ion density modulations whereas the predictions based on the ring model reproduce the 
simulation results accurately [22]. 

If the wake is excited by a self-modulated particle beam then it will grow in £ back 
through the bunch. If the S-M is constant then the wake will grow secularly in £. Under 
these conditions E rwake 3> -E r b ea m so the ponderomotive force can be obtained from Eq. 
We compared the predictions for rii using the full value for (E r ) in Eq. (JSJ) and the first 



term in the expansion e(E r ) ~ — e 2 /(4m e o; 2 )V r -E 2 for o z = 100c/co p , = 0.01 — 0.04n 
and oy <C 1. We assumed that E r grew secularly with a resonant driving term of E r f,. We 
find that keeping only the first term in the expansion underestimates the ion fluctuations, 
8 = (rii — n )/no in a wide range by 2-45%. 

As a guide, one can estimate the position within the beam for which the ion compression 
on-axis is smaller than an accepted value for 8. Assuming secular growth for E r , this occurs 
when £ cr it/az < 1 where £ cr i t is the required distance for 8 ~ 1 at r = 0. In self-modulated 
regimes the onset of the ion motion then occurs when: 

= f m,c 2 \ 1/2 f ^m^ _ 3eayE b r b \ 

a z \m e Z*lu*J V eVE b 8nm e c 2 1 r) ) 1 1 

Equation ^ shows that the ion motion can be minimized in the presence of heavier plasma 
ions (i.e. lower rrii/Z). The leading order term coincides with the fluid theory result, and the 
remaining terms are corrections associated with the generalized ponderomotive force. For 
the expected initial parameters of the PDPWFA, where n = 10 14 cm~ 3 , = 10 12 cm -3 , 




FIG. 1: Comparison between the ion density profile (An = n,j — no) analytical predictions (dashed 
lines) with 3D (a-b) and 2D (c-f) OSIRIS PIC simulations (solid lines). (a,c,e) 3D on-axis plasma 
proton density profile. (b,d,f) corresponding transverse density profiles at the £ of the vertical 
dashed lines. 

and a z ~ 12 cm, Eq. ^ shows that the ion motion is not important since S, lon /cr z ^> 1. 
However, for self-modulated beams £ C rit/c"z > 1, and the plasma ion motion can become 
important. To lowest order, the onset of ion motion beyond £ = a z can be reached by using 
more massive ions such that [mj/m e ] SM > 0.08<x 4 [cm]n|[10 12 cm~ 3 ]/<5. For a 5 = 0.05, a 
mass ratio of [mj/m c ] SM > 40 x 10 3 is required. This indicates that singly ionized Argon 
plasmas (or heavier) could be used to avoid the deleterious effects of the ion dynamics. 

We have performed fully self-consistent simulations of the PDPWFA using 2D cylin- 
drically symmetric geometry. The simulations used a moving window with dimensions 
680 x 8(c/ujp) 2 , divided into 13600 x 320 cells, with 2x2 electrons and ions per cell. 
Plasmas with mi/m e = 1836 (H + ), mi/m e = 73440 (Ar + ), and m 8 = oo (immo- 
bile ions) were considered with n = 10 14 cm -3 . An SPS-LHC CERN proton bunch 
was initialized with an energy of 450 GeV and with a half-cut density profile given by 
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Ub = fibo 1 + cos yy-n/2{z — z )/a z j exp [— r 2 /(2cr^)] for < z < z = er z V2vr, with 
nbo/n = 0.0152, o z = 225.6 c/uj p = 12 cm, and a r = 0.376 c/uj p = 200 jum. 

The simulations confirm that the motion of the plasma ions can be neglected until the 
bunch significantly self-modulates. Once the bunch self-modulates, the (H plasma) ion 
compression strongly modifies the wakefields. The plasma electron and ion density after the 
beam has been self-modulated are shown in Figs. [2^i-b. Ions have not moved significantly 
at the front of the proton bunch (£ < 250c/ u p ). However, for £ > 250c/ u p , the on-axis ion 
density grows such that Hi/no ^> 1. Ion density voids, where Hi/no — 0, are also formed 
near r ~ a r . The strong ion motion leads to electron trajectory crossing (Fig. [2^,). The 
electron flow becomes turbulent, the plasma reaches a quasi-neutral state, and the wakefield 
disappears for £ > 400 c/tu p (Fig. [3]). 

If a H plasma is used the SMI slows down or becomes suppressed at the back of the beam 
(£ > 250c/o;p) because the plasma becomes quasi-neutral and fields become smaller. This is 
why the proton bunch density modulations are less pronounced in Fig. [3Jd than in Fig. [3^i. 
Although the bunch is fully self-modulated in Fig. [3^, some of the protons are still inside the 
box as revealed by the structures at r ~ c/u p . The reason for the slowdown/suppression 
of the SMI is also seen in Fig. [3J3 showing that the SMI driving field E r — Bq vanishes 
at the back of the bunch. The amount of energy that can be transferred to accelerated 
particles is also smaller when ion motion occurs (inset of Fig. [3)3). The inset of Fig. |3]3 
also illustrates the slowdown of the SMI when plasma ion motion occurs. Similarly to the 
focusing force, and also because the plasma is quasi neutral, the accelerating fields drop 
abruptly for £ > 250c/ u p , i.e. when plasma ions move. This further reduces the proton 
bunch energy modulations, and the energy that can be transferred to accelerated particles. 
Little ion motion occurs when Ar + ions are used, and the SMI and -E a ccei develop as if ions 
are immobile. This is in excellent agreement with our theoretical results. These results 
strongly suggest that a future PDPWFA experiment could be performed with Ar plasmas. 
There are additional factors that may also change the growth of the instability including 
detuning due to strong energy chirps along the beam, or the presence of plasma density 
ramps [23] . In realistic conditions, however, these effects are not as important as the plasma 
ion motion. 

In conclusion, we showed that the plasma ion dynamics can strongly affect a future 
PDPWFA experiment. It causes the early saturation of the self-modulation instability, 
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FIG. 2: 2D cylindrically symmetric OSIRIS simulations of a PDPWFA at r = 11450/w p (6.1 m) 
in an Hydrogen plasma, (a) Plasma electron density. The inset shows sample electron trajectories 
colored according to the initial radius (b) Corresponding plasma ion density. The horizontal solid 
line represents the driving ion beam density profile. The arrow indicates the propagation direction. 
The vertical solid line is a lineout of rtj at ^ = 640c/w p . 

reduces the accelerating gradients, and hence limits the energy transfer from the driver 
to accelerated particles. The conditions to minimize the impact of the ion motion were 
identified. This work demonstrates that the deleterious effects associated with the plasma 
ion motion can be avoided by resorting to plasmas with higher ion charge to mass ratios. 

The authors acknowledge fruitful discussions with Prof. A. Pukhov. Work supported by 
the European Research Council (ERC-2010-AdG Grant 267841), by FCT (Portugal) through 
grants SFRH/BPD/71166/2010, PTDC/FIS/111720/2009, and CERN/FP/116388/2010; 
DE-FC02-07ER41500 and DE-FG02-92ER40727 US and NSF PHY-0904039 and PHY- 
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FIG. 3: OSIRIS simulation results of a PDPWFA at r = IU50/lo p (6.1 m). (a-b) show the 
proton bunch density, (c) Plasma focusing force profile with mobile plasma ions. The solid line is a 
lineout at r = 0.5c/ui p . The inset shows the evolution of the energy of E^. (d) on-axis accelerating 
gradients using mobile Ar + and H + ions. 



0936266. Simulations were done on the 1ST Cluster at 1ST, Jaguar supercomputer under 
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